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Executive  Summary 


As  marine  winds  strengthen,  wave  breaking  becomes  increasingly  widespread.  This  has 
important  consequences,  both  for  offshore  operations  and  geophysically,  yet  wave  breaking  is 
not  included  in  present  marine  forecasts.  This  project  addressed  the  challenge  of  incorporating 
wave  breaking  physics  into  coupled  sea  state  marine  weather  forecasting  models.  The  aim  was  to 
add  accurate  wave  breaking  and  sea  spray  forecasts  to  the  standard  suite  of  forecast  sea  state 
variables,  for  sea  state  conditions  ranging  from  light  to  very  severe. 

Our  approach  required  a  combination  of  observational  and  modeling  efforts.  We  synthesized 
results  from  our  field  data  analyses  into  accurate  parameterizations  for  wave  breaking 
occurrence  and  strength.  We  also  needed  to  modify  the  wind  input  source  term  for  compatibility 
with  the  other  source  terms  and  physical  constraints.  These  refinements  were  incorporated  in  the 
dissipation  and  wind  input  source  terms  in  the  spectral  wave  evolution  models  used  to  generate 
forecasts. 

Model  development  focused  on  wide  spectral  bandwidth  computations  of  the  directional  wave 
spectrum  and  its  spectral  tail  region.  We  used  an  ‘exact’  version  of  the  weakly  nonlinear  four- 
wave  interaction  source  term  in  the  radiative  transfer  evolution  equation  for  the  wave  spectrum. 
We  carefully  validated  cases  for  which  data  exists,  ensuring  that  the  modeled  wave  properties, 
wind  input  and  upper  ocean  dissipation  rates  are  all  consistent  with  observed  levels.  This  was  a 
key  focus  of  our  effort  since  matching  the  source  term  levels  over  a  wide  range  of  wind  speeds  is 
essential  for  reliable  forecasts  of  breaking  waves,  sea  spray  and  the  enhanced  air-sea  fluxes  they 
generate. 

The  model  output  results  were  evaluated  critically  in  systematic  tests  from  light  to  severe  wind 
speeds.  After  extensive  iteration  and  refinement,  the  model  framework  is  ready  for 
implementation  and  further  testing  in  operational  models.  Given  that  our  framework  used  the 
‘exact’  form  of  the  nonlinear  wave-wave  interaction,  we  anticipate  that  the  least  amount  of 
adaptation  and  best  results  will  be  obtained  if  the  new  TSA  (Two  Scale  Approximation)  source 
term  for  the  nonlinear  interactions  (Resio  and  Perrie,  2008)  is  implemented  operationally, 
replacing  the  standard  DIA  source  term. 

An  additional  goal  within  this  project  was  to  be  able  to  use  the  forecast  breaking  wave  properties 
to  compute  the  associated  sea  spray  flux  from  the  modeled  breaking  wave  properties  and  the 
wind  speed  rather  than  from  the  wind  speed  alone.  To  accomplish  this  goal,  a  parallel  laboratory 
study  (SPANDEX  II)  was  undertaken  to  parameterize  sea  spray  flux  dependences  on  breaking 
wave  properties. 

Overall,  when  this  project  commenced  in  FY10,  existing  knowledge  on  wave  breaking  was 
limited,  and  it  was  necessary  to  advance  knowledge  of  the  fundamental  physics  of  wave 
breaking  in  order  to  be  able  to  incorporate  this  process  reliably  into  forecast  models.  Very 
significant  progress  on  the  scientific  fundamentals  was  achieved  during  this  project  in  parallel 
with  our  modeling  effort.  The  model  framework  that  has  been  developed  for  forecasting  the 
breaking  wave  properties  is  ready  for  operational  implementation.  While  detailed  knowledge  of 
the  physics  of  very  short  breaking  waves  remains  incomplete,  modeling  of  these  scales  is 
feasible.  The  capability  of  including  wave  breaking  properties,  especially  those  of  the  dominant 
waves,  as  part  of  operational  sea  state  forecasts  is  now  available.  This  will  increase  both  the 
accuracy  and  utility  of  the  next  generation  of  operational  coupled  sea  state  and  ocean  weather 
forecasting  models,  particularly  in  severe  to  extreme  sea  states. 
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1.  BACKGROUND  /  PROBLEM  STATEMENT 

In  2010,  the  National  Oceanographic  Partnership  Program  (NOPP)  initiated  a  five-year  research 
project  entitled  ‘Improving  Wind  Wave  Predictions:  Global  to  Regional  Scales’.  This  project 
was  focused  on  improving  operational  wind  wave  modeling,  by  transitioning  new  science  into 
such  models,  and  by  developing  new  physics  parameterizations  for  such  models.  A  detailed 
strategic  overview  appears  in  the  paper  by  Tolman  et  al.  (2013)  which  describes  the  general 
goals  of  the  project,  the  science  and  operations  gaps  it  attempts  to  bridge  and  the  data  sets  and 
validation  techniques  for  verifying  the  updated  operational  models.  This  report  overviews  our 
contribution  to  this  NOPP  research  initiative. 


1.1  Wave  breaking  in  coupled  air-sea  interaction  modeling 

Under  strong  forcing  conditions,  breaking  waves  are  a  conspicuous  feature  of  the  wind-driven 
sea  surface,  appearing  as  whitecaps.  The  impact  forcing  of  large  breaking  waves  provides  the 
greatest  safety  challenge  to  shipping  and  coastal  structures.  Through  their  overturning  of  the  air- 
sea  interface,  breaking  waves  profoundly  influences  the  dynamics  and  thermodynamics  of  the 
boundary  layers  of  the  upper  ocean  and  the  marine  atmosphere.  Consequences  of  breaking  in  the 
upper  ocean  surface  layer  include  greatly  enhanced  observed  turbulence  dissipation  rates  in  the 
near-surface  region  (e.g.  Terray  et  al.,  1996;  Gemmrich  and  Farmer,  2004;  Sutherland  and 
Melville,  2015).  In  the  atmospheric  marine  boundary  layer,  increased  wave  form  drag  results 
from  the  separated  air  flow  over  breaking  waves.  This  is  accompanied  by  augmented  sea  spray, 
bubble  production  and  acoustic  underwater  noise,  as  well  as  enhanced  microwave  backscatter 
and  emissivity.  These  numerous  and  diverse  aspects  of  wave  breaking  are  described  in  greater 
detail  in  review  articles  (e.g.  Banner  and  Peregrine,  1993;  Melville,  1996;  Sullivan  and 
McWilliams,  2010;  Perlin  et  al.,  2013). 


Given  the  important  role  played  in  air-sea  fluxes  by  wave  breaking,  we  set  out  to  include  this 
widespread  process  more  explicitly  in  the  next  generation  of  fully-coupled  forecast  models  for 
the  atmosphere-ocean  system.  The  aim  of  this  project  was  to  conduct  the  research  needed  to 
address  this  challenge,  as  no  such  modeling  capability  was  available.  Forecasts  of  breaking  wave 
properties  provide  valuable  practical  wave  climate  information  for  sea  keeping.  The  approach 
pursued  was  to  develop  more  realistic  parameterizations  for  wave  breaking  occurrence  and 
strength,  and  to  validate  them  against  observations  in  test-bed  models.  The  end  goal  was  to  be 
able  to  refine  these  improvements  ready  for  use  in  an  operational  sea  state  forecasting  system, 
such  as  WaveWatch  III.  There  are  also  enhanced  scalar  air-sea  fluxes  associated  with  wave 
breaking,  including  sea  spray,  bubbles  and  greenhouse  gases.  Such  information  is  needed  to 
underpin  future  operational  versions  of  coupled  atmosphere-wave-ocean  models. 


1.2  Wind-wave  modeling 


In  present  wind-wave  forecast  models,  breaking  waves  only  appear  implicitly  as  part  of  the 
spectral  dissipation  source  term.  Otherwise,  they  receive  no  direct  quantification.  These  models 
compute  the  evolution  of  the  directional  wave  height  spectrum  F(k,0;  x,y,t)  under  the  resultant  of 
the  spectral  source  terms,  according  to  the  radiative  transfer  equation  (here  assuming  deep  water 
and  no  currents): 


dF 


§ 


•VF 


tot 


dt 


(1) 
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where  F  =  F(k,0)  is  the  directional  wave  spectrum,  cg  is  the  group  velocity.  The  total  source 

term  Stot  =  Sin  +  Sni  +  Sds,  where  Sm  is  the  atmospheric  input  spectral  source  term,  Sni  is  the 
nonlinear  spectral  transfer  source  term  representing  nonlinear  wave-wave  interactions  and  Sds  is 
the  spectral  dissipation  rate,  primarily  due  to  wave  breaking. 

In  present-generation  forecast  models,  the  wind  input  term  was  synthesized  from  ocean, 
laboratory  measurements  and  idealized  mathematical  models.  The  nonlinear  transfer  term  was 
derived  from  weakly  nonlinear  spectral  interactions  for  an  ensemble  of  irrotational  gravity  wave 
Fourier  modes  (Hasselmann,  1962).  The  dissipation  term  was  formulated  to  close  the  problem 
and  was  based  on  a  decay  rate  formulation  controlled  by  an  integral  spectral  wave  steepness 
parameter. 

It  should  be  noted  that  the  source  terms  in  these  models  are  heavily  tuned  to  produce  wave 
height  spectral  estimates  that  agree  well  with  corresponding  observations.  This  does  not  validate 
the  relative  strengths  of  the  three  source  terms,  which  are  more  difficult  to  assess.  One  option  is 
to  compare  computed  and  measured  integrated  wind  input  and  dissipation  rates,  for  which  the 
integrated  nonlinear  transfer  term  vanishes.  Unfortunately,  neither  the  integrated  wind  input  or 
dissipation  rate  is  computed  reliably  in  the  models,  nor  are  they  easily  measured.  Hence  this 
approach  has  not  been  actively  pursued  by  investigators.  However,  the  recent  availability  of 
systematic  breaking  wave  data  provides  a  more  robust  basis  for  discriminating  between 
alternative  wind  input  and  dissipation  rate  source  terms,  and  was  exploited  in  this  project.  Our 
model  validation  also  depends  on  detailed  comparison  of  forecast  wave  spectra  with  measured 
high  resolution  directional  ocean  wave  spectra,  which  have  become  available  in  recent  years 
(e.g.  Romero  and  Melville,  2010). 

1.3  Breaking  wave  physics  and  observations 

This  investigation  both  motivated  and  generated  several  significant  advances  in  breaking  wave 
observations  and  analysis.  Prior  to  this  project,  available  wave  breaking  field  data  showed  a  very 
large  scatter  when  plotted  against  standard  sea  state  variables  such  as  wind  speed  or  wave  age 
(e.g.  see  Figs.  7,  8  in  Gemmrich  and  Farmer,  1999). 

This  scatter  has  seen  significant  reduction  over  the  past  two  decades.  The  observational  study  by 
Banner  et  al.  (2000)  reported  a  strong  correlation  of  breaking  probability  of  dominant  wind 
waves  with  their  significant  steepness  above  a  threshold  level.  However,  an  extension  of  this 
class  of  result  to  the  shorter  waves  in  the  spectrum  was  needed.  Using  the  1999  Pacific  Ocean 
storm  waves  data  [Gemmrich  and  Farmer  (IOS,  Canada)],  Banner  et  al.,  2002)  found  a  common 
breaking  threshold  applicable  across  different  wave  scales.  This  threshold,  based  on  the  wave 
saturation  spectrum,  contributed  a  key  element  of  our  spectral  dissipation  rate  source  term. 

A  suitable  framework  was  needed  to  forecast  breaking  wave  properties,  which  was  provided  by 
the  spectral  density  of  breaking  crest  length/unit  sea  surface  area  (A(c))  introduced  by  Phillips 
(1985).  Measurements  of  A(c)  are  now  available  for  a  representative  range  of  wind  speeds  and 
wave  ages  through  remotely  sensed  imagery  (Kleiss  and  Melville,  2010;  Gemmrich  et  al.,  2013; 
Schwendeman  et  al.,  2014;  Sutherland  and  Melville,  2015).  These  results  have  been  of  primary 
importance  in  calibrating  the  source  terms  in  this  study.  A  detailed  description  of  A(c),  including 
available  ocean  measurements,  is  given  in  section  2.1.3.  While  A(c)  is  a  key  validation  variable, 
it  is  kinematic  in  nature.  A  breaking  strength  parameter  is  still  needed  in  conjunction  with  A(c), 
which  is  also  described  in  detail  in  section  2.1.3. 
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1.4.  Wind  stress  and  upper  ocean  dissipation  rate 

For  accurate  coupled  modeling,  the  wind  stress  and  upper  ocean  dissipation  rate  generated  by  the 
wave  model  need  to  be  consistent  with  observed  levels  and  their  trends  with  wind  speed  and 
wave  age.  The  recent  papers  by  Edson  et  al.  (2013)  and  Yang  et  al.  (2014)  for  ocean  wind  stress 
and  Sutherland  and  Melville  (2015)  for  the  upper  ocean  dissipation  rate,  provide  excellent 
observational  data  for  validating  the  compatibility  of  wave  forecast  model  capabilities. 


2.  TECHNICAL  APPROACH 

We  developed  a  modeling  framework  and  performed  numerical  model  calculations  for  wind- 
wave  spectral  evolution.  The  output  was  validated  against  a  set  of  benchmarks  based  on  field 
measurements.  The  source  terms  were  refined  iteratively  to  improve  the  model  performance. 
Due  to  the  complex  nature  of  the  ‘exact’  nonlinear  wave  interaction  source  term,  the  calculations 
typically  needed  24-48  hours  of  CPU  time  to  process  the  evolution  to  full  development.  Also, 
this  project  required  a  strong  synergy  between  our  modeling  effort  and  basic  knowledge  of  the 
physics  of  breaking,  including  the  analysis  of  field  observations.  This  was  an  essential  element 
that  underpinned  our  progress  throughout  this  project. 

2.1  Breaking  wave  observational  inputs 

Fundamental  observational  contributions  on  breaking  waves  generated  within  this  effort  were 
introduced  briefly  in  section  1.3,  and  are  highlighted  in  section  3.1.  Detailed  results  are 
described  in  sections  4.4  and  4.8. 

These  results  materially  advanced  the  development  of  our  goal  of  including  breaking  wave 
physics  in  our  model  framework.  Key  elements  were:  (i)  the  formulation  and  refinement  of  the 
spectral  dissipation  rate  source  term  Sds;  (ii)  the  extraction  of  spectral  breaking  wave  crest  length 
and  strength;  (iii)  incorporating  the  strong  enhancement  of  the  local  wind  momentum  flux  to  the 
waves  over  breaking  waves  into  the  wind  input  parameterization  term  (Banner,  1990a;  Kukulka 
and  Hara,  2008a, b). 

2.1.1  Breaking  onset  dependence  on  spectral  saturation 

A  collaborative  study  (Banner  et  al.,  2002)  undertaken  with  Gemmrich  and  Farmer  at  IOS 
Canada  involved  a  novel  scale  analysis  of  breaking  waves.  This  was  motivated  by  the  results  of 
the  model  study  of  wave  breaking  onset  by  Song  and  Banner  (2002),  who  identified  the  wave 
energy  convergence  rate  and  geometrical  steepening  within  nonlinear  wave  groups  as  key 
aspects  in  wave  breaking  onset.  The  directionally- normalized  spectral  saturation  <r(/)  was 
chosen  as  a  surrogate  for  local  wave  nonlinearity.  It  is  given  by 

cx(f)  =(27r)4f5G(f)/2g2/0(f)  (2) 

where  f  is  the  wave  frequency,  G(f)  is  the  wave  energy  spectrum  and  0(f)  is  the  directional 
spreading  width  spectrum. 

Banner  et  al.  (2002)  showed  that  a  common  breaking  threshold  based  on  a  spectral  wave 
steepness  parameter  collapsed  the  breaking  probability  data  over  an  appreciable  range  of  wave 
scales.  This  parameter  underpins  our  formulation  of  the  dissipation  source  term  and  our 
procedure  for  extracting  breaking  wave  properties.  More  specifically,  the  component  of  our  Sds 
formulation  associated  with  active  breaking  associated  with  a  given  wave  scale  bandwidth  was 
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built  around  the  observed  normalized  saturation  threshold.  Suitable  spectral  smoothing  was 
applied  to  the  measurements  of  normalized  spectral  saturation  <r . 

2.1.2  Spectral  dissipation  rate  source  term 

We  used  a  refined  version  of  the  saturation  threshold  form  of  the  dissipation  rate  source  term  Sds 
proposed  by  Alves  and  Banner  (2003),  incorporating  the  observed  breaking  saturation  threshold 
reported  by  Banner  et  al.  (2002).  Also,  various  refinements  were  introduced  to  shape  its  spectral 
distribution,  in  order  to  provide  a  much  better  match  to  the  wind  input  source  function  Sin  at 
higher  wavenumbers,  and  to  recently  published  observations  of  the  spectral  density  of  breaking 
wave  crest  length/unit  area  (Kleiss  and  Melville,  2010;  Gemmrich  et  al.,  2013;  Sutherland  and 
Melville,  2013). 

Sds  was  composed  of  a  local  breaking  term  S^'  plus  a  parametric  background  dissipation  rate 
Snd°nloc  quantifying  the  effects  of  longer  breaking  waves  attenuating  shorter  waves  (e.g.  Banner  et 
al.,  1989)  and  ambient  surface  layer  turbulence. 

The  form  of  Sdsused  in  this  study  is  shown  in  (3)  below.  It  uses  a  power  law  function  of  the 

normalized  spectral  saturation  ratio  to  reflect  the  observed  breaking  threshold  behaviour,  and 
refines  the  form  proposed  by  Alves  &  Banner  (2003).  The  form  used  here  is 

Sds(k ,9)  =  Sld°sc(k,6)  +  Srloc(k,0)  =C[(cr  -  crT)lcrT)a  +  £rJ(cr/crJb  «F(k,0)  (3) 

where  crand  a  are  respectively  the  saturation  and  saturation  normalized  by  the  mean  directional 
spreading  width,  crT  is  threshold  normalized  saturation  and  crm  is  the  saturation  at  km,  the  mean 
wavenumber  at  the  transition  from  the  peak  enhancement  region  to  the  spectral  tail. 

In  (3),  the  term  involving  the  normalized  saturation  threshold  ratio  (left  of  the  +  sign)  controls 
the  local  dissipation  rate.  The  breaking  threshold  switching  exponent  a  was  taken  as  2  and  b  was 
taken  as  2,  based  on  matching  to  the  high  wavenumber  forms  of  the  wind  input  source  term 
Sin(k),  as  discussed  below  in  section  3.2.1. 

The  coefficient  C  was  chosen  to  provide  the  optimal  match  to  observed  duration  evolution  data 
of  the  spectral  peak  energy  and  peak  frequency  (e.g.  Young,  1999).  It  is  evident  that  C  needs  to 
increase  with  the  wind  speed,  since  the  integrated  dissipation  rate  must  balance  the  integrated 
wind  input  for  developed  seas.  To  achieve  this,  we  used  a  function  of  the  mean  square  slope 
[mss]  of  the  wave  field  that  is  wind  speed  dependent  and  has  a  suitably  broad  dynamic  range. 

The  term  ^es  in  (3)  is  a  nonlocal  background  dissipation  coefficient  and  its  form  is  not  known. 

We  assumed  that  it  scales  with  the  mean  steepness  parameter  Ekp,  where  E  is  the  total  wave 

energy  and  kp  is  the  spectral  peak  wavenumber.  An  associated  proportionality  constant  of  0.125 
was  determined  through  model  runs  over  the  wind  speed  range  5-60  m/s,  for  which  its 
performance  was  found  to  be  effective.  This  parameterization  provides  attenuation  rates 
consistent  with  observed  decay  rates  of  swell  leaving  storm  areas  (Ardhuin  et  al.,  2009).  The 
measured  water-side  energy  dissipation  rate  s  in  the  wave  boundary  layer  is  the  sum  of  the  two 

contributions  Sldsc  (c)  +  Snd°nloc(c) ,  integrated  over  all  scales  and  over  the  depth  of  the  wave 
boundary  layer. 
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2.1.3  Spectral  density  of  breaking  crest  length  and  breaking  strength  parameter 
Spectral  density  of  breaking  crest  length/unit  area  A(c) 

As  introduced  in  section  1.3,  the  spectral  measure  of  breaking  wave  crest  length,  A(c),  has  the 
property  that  A(c)  dc  gives  the  crest  length/unit  surface  area  of  breaking  crests  traveling  with 
velocities  in  (c,  c+dc).  In  this  study,  A(c)  was  one  of  the  primary  breaking  forecast  parameters 
produced  from  our  modeling  effort  and  was  central  to  validating  our  source  terms  and  forecast 
breaking  wave  properties.  A  number  of  allied  air-sea  fluxes  are  based  on  the  sea  state  and  wind 
field  rather  than  on  the  wind  field  alone,  including  breaking  wave  enhancements  to  the  wind 
stress  and  input  to  the  waves,  as  well  as  sea  spray.  These  contributions  can  be  quantified  through 
modeled  breaking-related  properties. 

Spectral  breaking  strength  b(c) 

A  breaking  strength  measure  is  needed  to  quantify  the  wave  energy  loss  rate  through  breaking 
associated  with  a  given  A  distribution.  The  connection  between  these  two  distributions  was 
assumed  to  be  given  by  the  scalar  form  of  equation  (6.3)  in  Phillips  (1985),  which  we 
generalized  to  allow  the  breaking  strength  coefficient  to  have  a  spectral  dependence  [ b=b(c )]: 

Sdsc(c)  dc  =  b{c )  g~2c5A(c )  dc 

where  Sld°sc(c)  is  the  spectral  energy  loss  rate  from  wave  breaking,  as  defined  in  equation  (3). 

Underpinning  (4)  is  the  assumption  that  the  mean  wave  energy  dissipation  rate  at  scale  (c,  c+dc) 
is  due  primarily  to  wave  breaking  in  a  bandwidth  about  that  scale.  This  is  most  likely  valid 
around  the  spectral  peak,  however  A  levels  for  shorter  breakers  may  need  to  be  modified  due  to 
the  attenuation  of  short  wave  energy  by  longer  breaking  waves  (Banner  et  al.,  1989). 

It  is  noted  that  the  dependence  of  the  spectral  breaking  strength  coefficient  b(c)  is  not  well 
understood.  Initial  efforts  to  quantify  the  breaking  strength  were  based  on  the  Phillips  (1985) 
assumption  that  b  is  constant  for  all  breaking  wave  scales.  Measurements  were  made  of 
dissipation  rates  of  narrow-band  2D  laboratory  breaking  waves  in  focusing  packets.  These 
measurements  showed  a  strong  dependence  of  b  on  mean  packet  wave  steepness  or  energy 
convergence  rate  at  the  packet  envelope  maximum  (Banner  and  Peirson,  2007;  Drazen  et  al., 
2008;  Tian  et  al.,  2010).  While  measurements  of  b(c)  in  a  broad  wave  spectrum  are  not  presently 
feasible,  recent  parametric  forms  for  b(c)  have  been  proposed  that  are  based  on  a  spectral 
analogy  of  the  narrow-banded  measurement  findings  (Romero  et  al.,  2012;  Morison  and  Banner, 
2016).  These  are  formulated  as  a  function  of  the  normalized  characteristic  spectral  saturation. 
The  proposed  formulation  assumes  b(c)  ~  ( a  -  oT  )n,  where  a  and  !jt  are  the  normalized  spectral 
saturation  and  breaking  threshold  normalized  spectral  saturation,  and  n  is  an  exponent  to  be 
determined.  Further  details  and  results  based  on  this  approach  are  discussed  below  in  section  3.4. 

2.1.4  Effective  breaking  strength  beff 

Before  proceeding  to  review  our  model  results  against  of  our  proposed  suite  of  performance 
criteria,  we  note  that  the  spectral  breaking  strength  b(c)  has  certain  complexities,  as  mentioned 
above,  which  can  complicate  its  application.  Primarily,  b(c)  cannot  yet  be  measured  directly. 
The  only  relevant  measurements  presently  available  are  the  total  water-side  total  dissipation  rate 
as  a  function  of  depth,  hence  the  best  available  validation  method  is  indirect.  This  is  discussed  in 
detail  in  section  4.4  below. 

An  alternative  approach  is  to  use  a  mean  breaking  strength  coefficient,  beff,  defined  as 
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beff  =  j  SdT(c) dc  /g2  j  c5  *  A(c)  dc  (5) 

where  the  integrated  Sld°sc(c)  distribution  and  A(c)  are  either  obtained  from  wave  model  forecasts 

or  from  measurements.  The  variation  of  beff  with  wave  age  was  investigated  recently  by  us  in 
Zappa  et  al.  (2016),  based  on  measurements  and  spectral  wave  model  results.  This  is  reported  in 
section  4.4  below. 

2.2  Modeling  approach  and  objectives 

The  overall  objective  was  to  refine  our  model  framework  to  be  able  reproduce  the  field-observed 
benchmarks  A-E  stated  below  in  section  2.2.2  over  a  wide  range  of  wind  speeds.  Further,  this 
refined  modeling  framework  should  be  suitable  for  transitioning  to  operational  forecasting  in  the 
next  generation  of  coupled  atmosphere/sea  state/upper  ocean  circulation  models. 

2.2.1  Approach 
Wave  model 

We  developed  a  model  framework  with  a  spectral  dissipation  source  term  Sds  based  on  the 
threshold  nature  of  breaking.  This  term  was  assumed  to  depend  on  wave  parameters  and  is 
underpinned  by  the  saturation  threshold  behaviour  reported  by  Banner  et  al.  (2002).  The 
approach  was  based  on  treating  waves  in  different  directional  spectral  bands  as  nonlinear  wave 
groups.  It  belongs  to  the  class  of  nonlinear  forms  of  Sds  discussed  by  Donelan  and  Yuan  in  §11.4 
of  Komen  et  al.  (1994).  This  form  of  Sds  provides  a  method  for  calculating  A,  the  spectral 
density  of  breaking  crest  length/unit  area  (defined  above  in  section  2.1.3).  Details  of  this  form  of 
Sds  are  given  in  section  3  below.  In  addition  to  the  Sds  term,  we  needed  a  suitable  wind  input 
source  term,  having  committed  to  use  the  ‘exact’  form  of  nonlinear  spectral  transfer  term  S„i,  in 
preference  to  the  very  approximate  DIA  version.  Special  attention  was  given  to  ensuring  a  close 
match  between  the  total  energy  flux  to  the  waves  and  the  total  water-side  dissipation  rate  in  the 
wave  boundary  layer. 

Wind  forcing  model 

Our  strategy  in  this  project  was  to  utilize  an  accurate  wind  input  source  term  to  drive  the  model, 
with  an  appropriate  level  of  coupling  so  that  future  transition  to  a  fully-coupled  atmospheric 
boundary  layer  model  would  be  straightforward.  This  required  very  considerable  effort  and 
proved  to  be  very  challenging,  given  the  large  dynamic  range  of  growth  rates  in  existing  wind 
input  source  term  formulations.  In  fact,  this  project  was  a  test-bed  for  wind  input  source  terms, 
as  our  model  validation  procedure  placed  stringent  demands  on  this  source  term  for 
compatibility  with  the  other  source  terms  and  matching  to  observations,  as  outlined  in  section 
2.2.2.  Details  of  our  evaluation  procedure  and  the  final  form  of  our  wind  input  source  term 
implementation  are  summarized  in  3.2.1  below. 

2.2.2  Modeling  benchmarks 

A.  evolution  of  non-dimensional  mean  wave  energy  E  and  spectral  peak  frequency  fp. 

B.  spectral  tail  properties:  directional  lobing;  mean  directional  spreading  with  k/kp;  spectral 
saturation;  level  and  exponent  of  ID  transect  k-spectrum  (k  is  the  wavenumber,  kp  is  the 
spectral  peak  wavenumber);  frequency  spectral  transition  from  f4  to  f5. 
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C.  validating  the  relative  size  of  the  computed  wave-induced  stress  levels  (non-breaking  and 
breaking)  relative  to  overall  wind  stress  driving  the  model. 

D.  prediction  of  breaking  crest  length/unit  area  spectral  density  A  and  spectral  breaking 
strength  b(c)  of  the  wind  wave  spectrum  for  different  wave  ages. 

E.  ensuring  that  the  modeled  integrated  water-side  dissipation  rates  match  observed  levels,  and 
are  consistent  with  the  integrated  energy  flux  from  the  wind  to  the  waves. 

2.2.3  Spectral  bandwidth  of  the  calculations 

We  ran  broad  spectral  bandwidth  computations  of  the  evolution  of  the  directional  wave  spectrum 
and  its  tail  region  using  an  ‘exact’  version  of  the  nonlinear  wave-wave  interaction  source  term 
Sds  in  the  radiative  transfer  equation  for  the  wave  field.  Our  Sds  needed  to  be  uncompromised  by 
approximations  in  Sni.  We  used  the  ‘exact’  form  of  Sni  to  avoid  the  known  errors  associated  with 
‘discrete  interaction  approximation’  implementations  in  use  operationally. 

The  consequence  of  using  an  ‘exact’  form  of  S„i  was  very  long  computational  times  throughout 
the  project,  due  to  the  millions  of  interactions  that  need  to  be  calculated  at  each  space  or  time 
step  in  the  evolution.  This  effect  became  more  pronounced  at  higher  wind  speeds,  where  the 
time  or  space  step  needed  to  be  shortened  for  stability.  At  hurricane  wind  speeds,  model  runs  to 
well-developed  seas  required  several  days  to  finish.  Before  proceeding  with  the  calculations,  we 
invested  very  considerable  effort  in  evaluating  various  versions  of  the  ‘exact’  Sni  code  and 
propagation/stepping  schemes  available  to  us.  This  was  to  ensure  accuracy  and  minimize 
computational  instabilities  that  can  develop  at  higher  wavenumbers. 

For  our  standard  computations,  we  concentrated  on  duration- limited  growth  cases,  which  were 
the  least  prone  to  spurious  numerical  instability.  There  are  established  non-dimensional  duration- 
limited  growth  correlations  for  wave  energy  and  peak  frequency  based  on  field  observations 
(Young,  1999).  These  can  also  be  inferred  from  fetch- limited  observations  by  the  methodology 
discussed  in  Hwang  and  Wang  (2004).  In  our  calculations,  a  variable  logarithmic  grid  was 
developed  and  implemented  to  accommodate  the  evolution  from  very  young  to  very  old  wind 
seas  over  the  wide  range  of  wind  speeds  investigated. 

3.  TECHNICAL  ACTIVITIES  PERFORMED 
3.1  Summary  of  observational  data  analysis  effort 

(i)  we  prepared  a  compilation  of  all  published  breaking  probabilities  of  dominant  (spectral  peak) 
waves  against  normalized  spectral  saturation.  This  is  used  to  quantify  the  additional  breaking 
wave  stress  component  in  our  wind  input  source  function  and  to  provide  a  validation  check  of 
our  forecast  breaking  wave  spectral  crest  density  distributions. 

(ii)  Zappa,  Banner,  Morison  and  Brumer  (2016)  synthesized  from  the  available  open  ocean  data 
sets  a  compact  form  of  the  dependence  of  the  observed  mean  breaking  strength  coefficient  beff 
against  the  wave  age  based  on  data.  This  was  also  used  to  validate  our  wave  model  forecasts  of 
breaking  wave  properties  (section  4.4) 

(iii)  Fairall,  Zappa,  Banner,  Morison  and  Peirson  (2012)  conducted  a  follow-on  laboratory  study 
(SPANDEX  II)  which  significantly  refined  the  breaking  wave  measurements  gathered  in 
SPANDEX  I  (Fairall  et  al.,  2009),  which  was  the  first  study  seeking  to  relate  sea  spray  flux 
directly  to  breaking  wave  surface  disturbance  properties,  rather  than  wind  speed.  These  results 
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confirmed  a  strong  threshold  behavior,  with  a  quasilinear  increase  in  spray  flux  with  breaking 
energy  intensity  and  near-surface  wind  speed.  This  was  used  in  our  unified  sea  spray  flux 
forecast  model  formulation  described  in  4.8,  with  further  details  of  SPANDEX  in  section  4.1 1. 

(iv)  modeling  the  evolution  of  the  ocean  wind  wave  spectrum  depends  on  the  source  terms  and 
spectral  bandwidth  used  to  compute  the  governing  transport  equation.  Through  detailed 
measurement,  Sutherland  and  Melville  (2015)  investigated  the  relative  contributions  of  breaking 
wave  scales,  from  large-scale  whitecaps  to  micro-breakers,  to  the  total  dissipation  rate  in  the 
wave  boundary  layer.  They  concluded  micro -breakers  and  very  small  whitecaps  contributed  a 
“large  fraction  of  dissipation”  and  “that  a  large  fraction  of  wave  energy  was  dissipated  by  these 
small  breaking  waves.”  In  Banner  and  Morison  (2016),  we  reanalyzed  their  data  and  concluded 
that  for  young  seas,  micro -breakers  and  small  whitecaps  contribute  only  a  small  fraction  of  the 
total  breaking  wave  dissipation  rate.  For  older  seas,  micro -breakers  and  small  whitecaps 
contribute  a  large  fraction  of  the  breaking  wave  dissipation  rate,  but  this  is  only  a  small  fraction 
of  the  total  dissipation  rate,  which  is  now  dominated  by  the  non-breaking  component.  Hence,  for 
all  the  wave  ages  observed,  micro -breakers  and  small  whitecaps  do  not  appear  to  make  a 
significant  contribution  to  the  total  dissipation  rate  in  the  wave  boundary  layer. 

3.2  Model  development  and  computations 
3.2.1  Details  of  Source  Term  Development 
A.  Spectral  wind  input  source  term  Si„ 

We  used  a  form  of  the  wind  input  source  term  SM1  with  a  plausible  distribution  of  wind  input  to 
both  fast  and  slower  moving  waves  that  provides  accurate  energy  input  levels  from  weak  to 
hurricane  forcing  conditions,  consistent  with  observational  data.  We  also  investigated  the 
sensitivity  to  the  directionality  of  the  wind.  We  discovered  that  many  forms  for  Sin  are  non- 
optimal  in  our  modeling  framework  if  reliable  forecasts  are  sought  for  wave  breaking  properties. 

The  magnitude  and  spectral  composition  of  the  wind  input  source  term  are  imprecisely 
known,  despite  very  considerable  observational  and  theoretical  study  over  the  past  few  decades. 
In  BM10,  we  investigated  a  number  of  proposed  S,„  formulations:  Hsiao  &  Shemdin  (1983); 
Janssen  (1991);  Donelan  et  al.  (2006).  From  our  detailed  testing,  we  found  that  the  spectral 
distribution  of  the  growth  rate  was  crucial  to  the  successful  modeling  of  the  breaking  properties 
of  the  spectral  peak  waves.  With  too  low  a  level  of  input  to  the  dominant  waves,  there  was 
insufficient  excess  energy  flux  to  these  waves  to  account  for  the  observed  breaking  levels.  This 
key  factor  limited  out  attention  to  the  Janssen  (1991)  form  for  Sin,  which  was  modified  as 
discussed  in  detail  in  Banner  and  Morison  (2010)  through  a  tail- sheltering  modification. 

Conceptually,  the  wave  stress  driving  the  shorter  waves  is  reduced  from  the  total  wind  stress  by 
progressively  subtracting  the  wave  stress  associated  with  the  longer  waves.  This  approach 
allowed  fine  tuning  the  integrated  wind  input  energy  flux  to  closely  balance  the  integrated 
energy  loss  rate  due  to  breaking.  This  is  an  important  validation  check  for  the  modeling,  and  we 
found  that  our  sheltering  algorithm  provided  wind  stress  estimates  that  agreed  very  closely  with 
the  observed  levels  as  the  wind  sea  aged.  The  modified  Janssen  (1991)  wind  input  growth  rate 
was  used  in  S1T1  for  all  our  calculations. 

We  implemented  a  modified  form  of  the  Janssen  (1991)  wind  input  growth  rate  formulation, 
with  a  ‘sheltering’  term  designed  to  reduce  the  growth  rates  to  the  slower  moving,  shorter  wave 
components.  This  was  done  to  maintain  the  overall  aerodynamic  drag  coefficient  at  levels 
consistent  with  observed  levels. 
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Sheltering  strategy 


The  total  aerodynamic  drag  is  the  sum  of  the  wave  drag,  the  breaking  wave  drag  and  the 
tangential  stress,  i.e. 


Ttot  ~Tw+  Tbw  +  Ttm g 


(6) 


The  wave  stress  is  expressed  as: 

tw  =  jj  [/  F(k,  6 )  /  c]  kdkdO ,  (7) 

where  the  growth  rate  y  is  given  below. 

The  breaking  wave  stress  is  expressed  as: 

Tbw  =  J TiJk)dk  =  j Prfor(k)  *  ^>(k)  *  Tv(k )dk .  (8) 


Here,  Prbr(k)  is  the  breaking  probability  at  scale  k  and  Fr(k)  is  the  ratio  of  actual  crests  to 

Fourier  modes  at  scale  k.  For  details  on  how  Prbr  (k)  and  Fr(k)  are  defined  and  parameterized 
based  on  observations,  see  Banner  and  Morison  (2010). 

The  tangential  stress  parameterization  is  expressed  as: 

Ttang  =  A  pair  U(0  with  A  =  max[l0-5, 1/ (tx  U^,)] ,  where  tx  =  325,  t2  =  0.9. 

This  was  based  on  the  tangential  drag  coefficient  behaviour  reported  by  Banner  and  Peirson 
(1998).  They  observed  that  the  tangential  drag  coefficient  was  a  decreasing  function  of  the  wind 
speed,  largely  independent  of  the  wave  age.  We  developed  a  parameterization  that  conformed  to 
their  data  trend  at  moderate  wind  speeds,  and  asymptoted  to  a  residual  level  of  0.0001  for 
hurricane  winds  of  50  m/s.  This  is  based  on  the  assumption  that  no  matter  how  strong  and 
widespread  the  air  flow  separation  becomes,  the  wind  always  maintains  some  residual  re¬ 
attachment  to  the  water  surface,  where  the  tangential  stress  will  be  non-zero. 


In  the  present  modeling,  the  overall  friction  velocity  is  given  byn*  = 


.  Also,  a  reduced 


friction  velocity  felt  by  the  nth  wavenumber  kn,  which  reflects  the  sheltering  of  the  short  waves 
by  the  longer  waves,  is  given  by 


red 


(U  =  jT,o,  ~ askXM0  +  hJO) 


for  1  <  n  <  N 


i=i 


(9) 


where  N  is  the  total  number  of  wavenumber  grid  points.  At  each  timestep,  the  model  recalculates 
the  total  and  local  u*  levels,  which  are  then  fed  back  into  the  wind  input  source  term.  This  key 
aspect  of  the  wind  input  source  term  is  described  in  the  following  section.  Here  axj,  is  a  constant 
that  reflects  the  level  of  sheltering.  This  parameter  allows  setting  the  sheltering  level  to  obtain  an 
optimal  match  between  the  wind  stress  generated  in  the  model  and  the  observed  level  as  a 
function  of  wind  speed  and  wave  age. 


Modified  Janssen  (1991)  wind  input  growth  rate  formulation 
Following  Janssen  (1991),  we  define: 

p  =  0 uJkc f  (gz0  / ul)exp(J xkI(u*  cos 61c)2),  (10) 

where  k=  0.4  is  the  Karman  constant,  Ji  =0.99,  c  is  the  phase  speed  and  //is  the  direction  of  the 
waves  relative  to  the  wind. 
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The  Miles  parameter  (Miles,  1957)  is  given  by: 

/?  =  /2//(ln(//))4/x'2,  for  jU<\ where  J2  =1.6.  (11) 

(3  -  0  for  jU>  1 . 

The  spectral  growth  rate  is  then  given  by: 

yik,  0)  =  £(3co  (u7d  (k)  /  cf  cos  0 ,  (12) 

where  s  is  the  ratio  of  air  to  water  densities,  and  the  corresponding  wind  input  source  term  is: 

S  m=yF(k,0).  (13) 

Fig.l  below  illustrates  the  typical  impact  of  the  sheltering  in  relation  to  the  Janssen91  growth 
rate  as  well  as  to  other  proposed  growth  rate  formulations  (Snyder  et  al.,  1981;  Hsiao  and 
Shemdin,  1983;  Plant,  1990),  for  reference.  Note  that  the  form  of  the  Janssen  (1991)  growth  rate 
parameterization  has  been  largely  followed,  with  the  Miles  coupling  parameter  based  on  the 

overall  u*/c,  but  input  to  the  shorter  waves  modified  by  using  u£ed  in  the  quadratic  forcing  term. 


inverse  wave  age  1 .0,  U1 0=1 2  m/s 


Figure  1.  This  logarithmic  plot  highlights  the  considerable  differences  between  the  spectral  growth 
rate  P  of  selected  proposed  forms  of  Si„.  The  differences  are  shown  for  maturing  seas  (Uio/cp~1.0). 
The  modified  Janssen91  curve  shows  the  extent  of  sheltering  introduced  for  the  slower  moving, 
shorter  waves  needed  to  bring  the  computed  wind  stress  into  agreement  with  observed  levels. 

B.  Spectral  nonlinear  transfer  rate  source  term  S„i 

The  version  of  ‘exact’  Sni  used  was  a  recent  update  (Don  Resio,  private  communication)  of 
Tracy  and  Resio  (1982)  with  ±180°  directional  coverage.  Our  decision  to  run  a  wide  spectral 
bandwidth  required  careful  testing  for  high  wavenumber  instabilities  and  was  expensive  in  both 
development  and  computational  time.  The  instabilities  were  minimized  by  using  short  timesteps, 
and  initial  dampening  of  rates  of  change.  Also,  energy  removal  at  high  wavenumber  end  of  the 
grid  achieved  by  imposing  a  k-4  tail  termination  on  the  computed  solution,  thereby  suppressing 
any  accumulation  at  the  end  of  grid  associated  with  Sni.  The  model  output  verified  that  the 
nonlinear  transfer  term  had  zero  net  integral  at  all  times. 
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C.  Spectral  dissipation  rate  source  term  Sas 

The  form  of  Sdsused  in  this  study  was  described  in  detail  in  section  2.1.2,  and  was  stated  as 
equation  (3): 

Sds(k,0)  =  C[(<r-CT T)/aTy  +£rJ(<j/am)hcoF(k,6>)  (14) 

The  power  law  function  of  the  normalized  spectral  saturation  ratio  reflects  the  observed  breaking 
threshold  behavior.  Where  the  normalized  spectral  saturation  falls  below  the  threshold,  the  local 
dissipation  no  longer  arises  from  wave  breaking.  This  proposed  form  (3)  for  Sds  refines  the 
integral  wave  steepness  threshold  used  in  the  quasilinear  form  of  Sds  presently  used  in  many 
operational  sea  state  forecast  models. 

3.3  Wind  and  wave  model  computations 

We  assumed  a  neutrally  stable  logarithmic  atmospheric  boundary  layer  structure  over  the  sea. 

U(z)  =  —  log(f)  (15) 

K  z° 

The  wind-wave  model  was  initiated  with  a  JONSWAP  spectrum  at  a  suitably  short  fetch, 
depending  on  the  wind  speed.  Changes  in  the  background  wind  field  due  to  the  waves  impacted 
the  wind  stress  that  forced  the  wave  model.  This  was  subject  to  numerical  control  during  the  spin 
up  period. 

3.3.1  Computations  of  duration-limited  wind  wave  evolution 

Computations  of  the  directional  wave  spectrum  were  made  for  the  spectral  bandwidth  covering 
gravity  waves  from  0.01  Hz  up  to  5  Hz,  using  the  source  terms  described  above.  We  focused  our 
development  initially  on  the  case  of  a  steady  forcing  wind  speed  of  Uio  =  12  m/s.  This 

corresponded  to  available  observations  where  data  on  winds,  waves  and  breaking  waves  were 

gathered  simultaneously  and  could  be  used  for  detailed  validation.  Once  this  case  had  been 
optimized,  model  evolution  runs  were  carried  out  for  wind  speeds  from  3-60  m/s,  and  ran 
stably  up  to  lOOm/s. 

Benchmark  D  (section  2.2.2)  proved  to  be  particular  strategic  significance.  This  involved  a 
comparison  of  forecast  and  observed  breaking  wave  properties  during  growing  seas  especially  at 
the  spectral  peak  where  the  proposed  relationship  (4)  is  most  likely  to  be  valid.  Not  only  does 
this  provide  a  tighter  constraint  on  the  form  of  the  spectral  dissipation  rate  source  term,  but  it  has 
the  benefit  of  reducing  the  uncertainty  in  the  form  of  Sm,  as  discussed  below. 

We  implemented  and  refined  our  modeling  strategy  for  extracting  the  relevant  wave  breaking 
parameters,  evaluating  the  wind  stress  components  and  coupling  the  wind  and  wave  models  so 
that  a  feedback  was  operative  in  which  the  computed  wind  stress  forced  the  wave  model 
interactively.  Concurrently,  we  investigated  breaking  wave  properties,  validating  the  results 
using  newly  available  breaking  wave  data. 

3.4  Extraction  of  breaking  wave  properties: 

3.4.1  Breaker  crest  length  per  unit  area  spectral  density  A(c)  and  breaking  strength  b(c) 

Banner  and  Morison  (2010)  described  in  detail  how  b  and  A  at  the  spectral  peak  can  be  extracted 
from  the  numerical  calculations  of  S'd°sc  (k)  and  normalized  spectral  saturation  a .  In  this  project, 

we  refined  this  capability  to  compute  b(c)  and  A(c)  over  the  full  spectral  scale  bandwidth 
resolved  in  the  model. 
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We  followed  an  approach  similar  to  Romero  et  al.  (2012)  which  mimics  spectrally  the  observed 
breaking  strength  dependence  on  wave  steepness  reported  in  laboratory  breaking  wave  studies  of 
Drazen  et  al.  (2008),  but  uses  our  normalised  spectral  saturation  to  mimic  this  in  the  spectrum. 
After  adjustment  to  obtain  the  equivalent  spectral  form  using  our  open  ocean  normalized 
saturation  breaking  threshold,  we  defined  the  breaking  strength  parameter  as: 

b(c)  =  H(5 - aT )  *  Abr  *(g°5-  ct°5)  (16) 

where  Abr  =  0.0025  is  a  breaking  strength  coefficient.  The  corresponding  breaking  crest  length 
distribution  then  follows  from  (5).  Hence 

A(c)=Si“(c)*g2/(b(c)*c!)  (17) 

Thus  A(c)  is  derived  from  our  modeled  S^'  (c)  and  (16)  for  b(c),  so  verifying  A(c)  validates 
S^(c)  and  b(c).  Note  that  At,r  is  determined  empirically  through  matching  modeled  and 
observed  B  effective  levels. 

This  post -processing  spectral  breaking  formulation  allows  the  breaking  strength  and  crest  length 
distributions  of  the  waves  to  be  forecast  from  the  output  of  a  spectral  wave  model.  However, 
there  are  potential  sensitivities  that  need  to  be  taken  into  account  as  they  can  affect  the  accuracy 
of  the  extracted  breaking  properties.  They  are  associated  with  the  rapid  variation  of  the 
normalized  spectral  saturation  <r  and  local  dissipation  rate  in  the  vicinity  of  the  spectral  peak, 
together  with  the  spectral  dependence  of  b(c). 

As  the  project  progressed,  we  were  able  to  significantly  refine  our  capability  for  calculating  the 
spectral  density  of  mean  breaking  wave  crest  length/unit  area  A(c)  to  the  point  where  close 
agreement  was  achieved  between  modeled  and  observed  levels  of  breaking  probability  for  the 
dominant  waves  for  both  developing  and  maturing  seas  (see  section  4.4  below).  Also,  the  recent 
availability  of  microscale  breaking  data  from  infrared  measurements  (Sutherland  and  Melville, 
2013,  2015),  provided  better  detection  and  resolution  of  the  very  short  breaking  waves  in  the 
spectral  wavenumber  tail  which  may  not  entrain  air.  Our  model  results  are  able  to  capture  their 
behavior  reasonably  well,  as  discussed  in  section  4.4  under  ‘Breaking  wave  forecasts’. 

3.4.2  Spectral  peak  wave  breaking  probability 

We  were  able  to  verify  the  model  forecasts  of  A(c)  independently  of  b(c)  by  comparing 
dominant  wave  (spectral  peak)  breaking  probabilities  calculated  from  the  forecast  A(c)  against 
observed  levels,  and  validate  our  best  fit  parameterisation  for  breaking  probability  against 
normalized  spectral  peak  saturation,  using: 


1-3  cp 

Prbr(c)  =£7  ” cA(c)dc/  | cll(c)dc 


0.7  c„ 


(18) 


where  cp  is  the  spectral  peak  phase  speed.  This  represents  the  ratio  of  breaking  crest  passage  rate 
to  the  passage  rate  of  all  crests  in  the  spectral  peak  region.  In  Appendix  A  in  Banner  and 
Morison  (2010),  it  is  shown  that  the  denominator  is  well-approximated  by  0.6  g/(27icp2).  From 
the  modeled  spectral  evolution,  the  spectral  peak  breaking  probability  can  be  computed  for  any 
wind  and  wave  age,  and  the  results  compared  with  the  ensemble  of  currently  available 
observations.  The  breaking  probability  correlation  against  normalized  saturation  for  the  data  we 
compiled,  as  well  as  our  model  parameterization  used  in  calculating  the  breaking  stress  is  shown 
in  Fig.2,  below.  Validation  of  model  results  for  breaking  properties  against  this  data  calculated 
using  equation  (18)  is  discussed  in  section  4.4  below. 
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Figure  2.  The  parametric  relationship  (solid  black  curve)  used  in  our  model  to  estimate  the 
additional  breaking  wave  stress  in  the  computed  wind  stress,  was  obtained  as  a  fit  to  the  measured 
data  (colored  circles). 

4.  WAVE  MODEL  EVOLUTION  RESULTS  AND  DISCUSSION 

We  began  by  validating  our  model  results  for  wave  energy  and  spectral  peak  frequency  against 
the  duration-limited  data  trend  curve  given  by  Young  (1999,  §5.3.4)  (benchmark  A  in  section 
2.2)  at  Uio=12  m/s.  The  model  diagnostics  for  Uio=12  m/s  were  very  encouraging,  so  we  carried 
out  model  runs  over  a  very  wide  range  of  wind  speeds  from  3-60  m/s  to  investigate  the 
robustness  of  our  model  performance.  It  is  important  to  note  that  the  model  settings  were  not 
modified  for  the  higher  wind  speed  cases  explored,  nor  were  any  limiters  used  to  constrain  the 
results.  Thus  benchmark  A  was  met  in  terms  of  closely  matching  the  evolution  of  the  standard 
growth  curves  over  the  entire  wind  speed  range  investigated  (3  to  60  m/s),  as  seen  in  Fig. 3. 

4.1  Benchmark  A  -  Evolution  of  total  wave  energy  and  peak  frequency  from  3-60  m/s 


Figure  3.  Evolution  of  non-dimensional  mean  wave  energy  (upper  panel)  and  spectral  peak 
frequency  (lower  panel)  against  non-dimensional  time  for  duration  limited  growth  under  3-60  m/s 
wind  forcing.  The  background  dashed  lines  are  the  trends  of  the  data  collated  by  Young  (1999). 
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4.2  Benchmark  B  -  Spectral  tail  properties 

The  results  for  benchmark  B  in  section  2.2  are  shown  in  Figs.  4  and  5.  These  spectral  measures 
have  been  checked  carefully  against  available  data  and  are  in  close  agreement. 


Figure  4.  One-dimensional  (ki)  wavenumber  transect  spectra  in  the  upwind-downwind  direction, 
for  wind  speeds  from  3  to  60  m/s  for  a  range  of  wave  age  conditions.  The  dashed  black  line  is  k'3. 


One-dimensional  transect  wavenumber  spectrum 

Computed  one-dimensional  (ki)  wavenumber  transect  spectra  in  the  upwind-downwind  direction 
are  shown  in  Fig.4  for  Uio  ranging  from  3  to  60  m/s  for  a  range  of  wave  age  conditions.  These 
spectra  show  a  close  correspondence  with  the  measured  kf3  tail  exponent  and  weak  wind  speed 
dependence  reported  by  several  studies  (e.g.  Romero  and  Melville,  2012). 

Directional  spreading  properties 

The  directional  spreading  of  the  spectrum  is  an  important  feature  that  has  been  measured  reliably 
in  recent  years,  and  is  an  important  feature  for  source  term  balance  validation.  Fig.  5  shows  the 
directional  spreading  at  different  k/kp  for  a  range  of  wind  speeds  (12  m/s.  24  m/s,  36  m/s  and  48 
m/s  and  two  wave  ages,  corresponding  to  young  seas  and  to  mature  seas. 

Characteristically  single-peaked  near  the  spectral  peak,  these  directional  spreading  distributions 
develop  a  distinctive  bi-modal  peak  towards  higher  k/kp  with  maxima  exceeding  up  to  +50 
degrees  off  the  primary  wave  direction,  and  lobe  ratios  up  to  1.2.  The  variability  with  respect  to 
k/kp,  wind  speed  and  wave  age  can  be  assessed  from  these  results  and  compared  with  available 
data.  The  comparison  is  generally  very  favourable,  but  maximum  measured  lobe  ratios  tend  to  be 
higher  than  our  modeled  levels.  The  Sni  term  is  at  the  heart  of  the  lobed  spreading  distribution, 
and  our  results  suggest  that  the  form  of  Sni  used  may  need  some  refinement  in  order  to  match  the 
lobe  ratio  data  more  closely.  However,  this  is  regarded  as  of  secondary  importance. 
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Figure  5.  The  eight  panels  show  directional  distribution  of  wave  energy  with  distance  from  the 
spectral  peak  for  12,  24,  36  and  48  m/s  forecasts.  The  left  panels  are  for  younger  seas,  and  the  right 
panels  are  for  older  seas. 

Mean  directional  spreading 

The  directional  spreading  angle  is  defined  as: 

0  =  J  6F(k,  6)kdkd0 /  J  F(k,  O)kdkd0 


(19) 
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Fig.  6  shows  the  modeled  variation  of  the  mean  directional  spreading  angle  with  distance  from 
the  spectral  peak  for  wind  speeds  3-60  m/s  for  young  seas,  mature  seas  and  for  a  full  range  of 
wave  ages.  It  also  shows  these  results  for  4-16  m/s  wind  speeds  and  wave  ages  against  data  from 
Romero  et  al.,  2010.  The  last  panel  shows  the  close  agreement  between  forecasts  and  data. 

(a)  (b) 


(c)  (d) 


Figure  6.  Modeled  variation  of  the  mean  directional  spreading  angle  with  distance  from  the 
spectral  peak,  (a-c)  for  wind  speeds  3-60  m/s  for  (a)  young  seas;  (b)  mature  seas  and  (c)  for  a  full 
range  of  wave  ages,  (d)  mean  spreading  angle  for  4-16  m/s  matching  wind  speed  and  wave  ages 
with  data  from  Romero  et  al.,  2010,  with  the  matching  data  plotted  in  heavier  black-dashed  lines. 

Normalized  spectral  saturation 

The  normalized  spectral  saturation  a  defined  in  section  2.1.1  is  used  in  the  formulations  of  the 
local  dissipation  and  breaking  strength  parameters.  Fig.7  below  shows  its  variation  with  k/kp  for 
the  wind  speed  ranging  from  3-60  m/s  for  two  wave  age  cases,  corresponding  to  young  and 
mature  wind  seas.  Also  shown  is  the  breaking  threshold  level  of  a  found  by  Banner  et  al.,  2002 
and  used  in  our  modeling. 

It  is  noted  that  this  parameter  is  akin  to  the  classical  ‘ak’  slope  steepness  parameter  that  serves  as 
a  measure  of  nonlinearity  in  wave  analysis.  While  a  geometric  parameter,  it  is  the  closest 
surrogate  that  has  been  proposed  to  date  for  parameterizing  the  nonlinear  dynamics/energetics  of 
wave  trains  in  a  spectral  context. 
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Figure  7.  Plots  of  the  behavior  of  the  normalized  saturation  against  k/kp  for  wind  speeds  3  to  60 
m/s.  The  red  dashed  lines  shows  the  threshold  used  in  the  local  dissipation  source  term.  Left  panel 
is  for  young  seas,  right  panel  for  older  seas. 


Frequency  spectrum  tail  transition 


Frequency  Energy  Spectra 


Frequency 


Figure  8.  Computed  frequency  spectra  at  different  stages  of  development  for  Uio=12  m/s.  Red  and 
blue  dashed  lines  have  exponents  -4  and  -5  respectively. 


The  frequency  spectra  computed  from  the  modeled  directional  wavenumber  spectra  show  the  tail 
slope  transitioning  from  f4  to  f5  at  about  3  times  the  spectral  peak  frequency.  This  behavior  is 
also  seen  in  many  wave  measurements  using  single-point  instruments  such  as  buoys  or  wave 
poles.  This  phenomenon  is  due  to  the  increasing  directional  spreading  with  distance  from  the 
spectral  peak,  as  first  explained  by  Banner  (1990b). 

4.3  Benchmark  C  -  Associated  wind  stress/drag  coefficient  properties 

Benchmark  C  (section  2.2.2)  is  concerned  with  the  relative  size  of  wave-induced  stress  level. 
The  sum  of  the  wave  stress  (non-breaking  and  additional  breaking  contributions)  and  the  viscous 
tangential  stress  equals  the  total  wind  stress.  The  most  recent  comprehensive  open  ocean  data  for 
measured  sea  surface  drag  coefficients  (Edson  et  al.,  2013)  covers  the  wind  speed  range  5-27 
m/s  and  a  subset  of  wave  age  conditions  of  the  computations. 
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Figure  9.  Behavior  of  the  sea  surface  drag  coefficient  with  wave  age  for  wind  speeds  from  3-60  m/s 
during  the  evolution  from  young  to  old  wind  sea  conditions.  The  superimposed  data  show 
measured  sea  surface  drag  coefficients  from  Edson  et  al.  (2013).  The  color  code  for  wind  speed  is 
seen  in  the  legend.  The  left  panel  has  ash  =0.95  and  the  right  panel  has  ash  =  0.85. 

The  model  results  shown  in  Fig. 9  agree  closely  with  these  observed  levels  and  indicate  likely 
levels  when  extrapolated  beyond  the  measured  ranges.  These  two  figures  illustrate  the  sensitivity 
to  the  size  of  the  sheltering  coefficient  ash-  It  is  seen  that  reducing  ash  from  0.95  to  0.85  improves 
the  agreement  between  modeled  and  observed  trends,  but  appreciably  increases  the  drag 
coefficient  for  very  high  wind  speeds.  As  there  is  no  drag  coefficient  data  available  yet  for  these 
extreme  winds  for  conventional  large  fetch  (non-hurricane)  conditions,  the  optimal  choice  of  ash 
remains  an  open  question. 

Normal  and  breaking  wave  stress  contributions  to  the  wind  stress 

As  discussed  above  in  section  3.2.1,  the  total  wind  stress  is  modeled  as  three  additive 
components:  the  wind  stress  associated  with  wind  momentum  flux  to  the  waves  (‘wave  stress’), 
a  wave  stress  enhancement  associated  with  separated  wind  flow  over  breaking  waves,  plus  a 
viscous  tangential  stress  component,  usually  small.  The  modeled  behavior  of  the  first  two 
contributions  is  shown  in  Fig.  10  below,  which  indicates  their  relative  importance  at  different 
stages  of  development  for  wind  speeds  ranging  from  3-60  m/s. 


Figure  10.  Variation  with  wave  age  of  the  normalized  wave  stress  (left  panel)  and  additional 
breaking  wave  stress  (right  panel)  for  the  wind  speed  ranging  from  3-60  m/s.  These  combine  with 
the  viscous  stress  to  obtain  the  total  momentum  flux  to  the  sea  surface  from  the  wind. 
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4.4.  Benchmark  D  -  Breaking  wave  forecasts 
Spectral  density  of  breaking  crest  length/unit  area 

A  major  aspect  of  our  model  validation  sought  to  reproduce  observed  breaking  wave  statistics. 

The  most  recent  comprehensive  open  ocean  measurements  of  A(c)  were  reported  by  Kleiss  and 
Melville  (2010),  Romero  et  al.  (2012)  and  Sutherland  and  Melville  (2015).  For  the  slower- 
moving  shorter  breakers,  previous  visible  video  imagery  data  showed  an  unexplained  strong 
attenuation  of  A  below  c=3  m/s.  However,  the  recent  resolution  and  sensitivity  advances  in 
infra-red  video  capability  for  imaging  breakers,  including  non-aerating  microscale  breakers.  This 
has  allowed  extending  the  slow-moving  breaker  range  significantly  (Sutherland  and  Melville, 
2013,  2015).  So  the  higher  resolution  A(c)  observations  more  closely  match  the  shape  and  level 
produced  by  the  model.  Further  enhancements  in  measurement  techniques  may  reduce  the 
difference  in  A(c)  for  speeds  below  1  m/s. 
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Figure  11.  Breaking  crest  length  per  unit  area  of  sea  surface  [A(k)]  versus  k/kP  for  young  seas  (left) 
for  wind  speeds  3  to  60  m/s  and  for  mature  seas  (right). 

The  A(k)  measurements  show  consistently  high  levels  of  A  in  the  tail,  and  an  enhanced  A 
around  the  spectral  peak  enhancement  region.  Also  as  the  seas  mature,  the  spectral  levels  of  A  in 
the  spectral  peak  region  reduce  significantly.  However,  to  date,  no  reliable  observations  of  A(k) 
are  available,  and  so  we  will  look  at  A(c)  and  b(c)  in  the  wave  speed  domain,  where  reliable 
A(c)  measurements  have  been  made. 


Figure  12.  Breaking  strength  b(c)  versus  wave  speed  c,  for  different  wave  ages  (Uio/cp)  for  12  m/s 
(left)  and  48  m/s  (right). 


The  breaking  strength  is  weakly  dependent  on  wind  speed,  but  significantly  dependent  on  wave 
age  and  wave  speed.  There  are  no  observations  for  b(c),  largely  because  it  is  difficult  to  measure 
the  speed  (or  wavenumber)  dependence  of  the  local  wave  dissipation.  We  can  validate  the 
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effective  breaking  strength  in  observations  using  the  integrated  dissipation  and  integrated  A 
distributions,  and  this  is  shown  in  Fig.  16  below. 

A  number  of  observational  data  sets  have  been  collected  for  A(c),  and  recently  this  has  been 
extended  into  the  shorter  waves,  using  infrared  video  by  Sutherland  and  Melville  (2015).  Fig. 13 
below  shows  the  forecasts  of  A(c)  from  the  model,  matching  wind  speed  and  wave  age  ranges 
with  the  data  from  Sutherland  and  Melville  (2015). 
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Figure  13.  Plots  of  breaking  crest  length  per  unit  area  (A(c))  versus  wave  speed  (c)  for  model  runs 
4-16  m/s  matching  the  wind  speed  and  wave  age  range  of  Sutherland  and  Melville  (2015).  The 
matching  observational  data  is  plotted  in  thick  black-color  dotted  lines,  with  the  colors  of  both  the 
observations  and  the  model  data  matching  the  wave  ranges  from  Sutherland  and  Melville  (2015). 


The  results  indicate  very  encouraging  agreement  between  observed  and  forecast  levels  during 
developing  sea  state  conditions.  For  ‘fully-developed’  seas,  the  model  produces  no  breaking  at 
the  spectral  peak.  The  observations  start  to  roll  away  from  the  c-6  slope  below  about  2m/s,  and 
the  roll  over  totally  below  1  m/s.  While  the  model  slope  similarly  decreases  below  2  m/s,  the 
model  A(c)  never  rolls  over.  However,  this  discrepancy  is  not  energetically  significant.  This 
difference  may  be  due  to  deficiencies  in  either  the  very  tail  of  the  model,  or  the  tail  of  the 
observations,  and  further  enhancements  in  observation  capability  or  in  our  understanding  of  the 
physics  of  short  breaking  waves  may  remove  this  minor  discrepancy  in  the  future.  Clearly, 
further  comparisons  with  data  are  now  needed  to  establish  the  robustness  of  the  modeling 
approach  proposed  here. 

While  further  field  observational  validation  is  required,  the  agreement  of  these  initial  results  is 
very  reassuring.  To  the  extent  possible,  this  phase  of  our  effort  met  our  modeling  benchmark  D. 


Breaking  probability  data  analysis  results 

Banner  et  al.  (2002)  reported  an  analysis  of  breaking  probability  with  scale  based  on  using  the 
directionally- normalized  spectral  saturation  cr(f)  =  (27i)4f  5G(f)/2g2/0(f),  where  f  is  the  wave 
frequency,  G(f)  is  the  wave  energy  spectrum  and  0(f)  is  the  directional  spreading  width 
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spectrum.  Their  results  revealed  a  remarkable  near-uniform  collapse  for  the  breaking  probability 
in  different  spectral  bands  at  and  above  the  spectral  peak  waves.  They  highlighted  a  common 
threshold  for  <r(f)of  approximately  4.5xl0'3  for  breaking  onset  and  established  order  into  the 
analysis  of  breaking  wave  measurements  for  the  first  time. 


Figure  14.  Breaking  probability  for  the  dominant  wind  sea,  calculated  using  equation  (18),  versus 
wave  wage  for  Uio  =  3  to  60  m/s.  The  color  coding  refers  to  the  normalised  saturation  (Nsat=(j ) 
shown  in  the  legend. 

As  described  in  detail  above  in  section  3.2.1,  our  present  study  used  a  formulation  for  Sds  for  the 
dissipation  rate  associated  with  wave  breaking  based  on  this  observed  normalized  saturation 
threshold.  As  part  of  the  wind  input  source  term,  there  is  a  component  associated  with  additional 
wave  form  drag  over  breaking  waves  associated  with  the  separated  flow  they  induce.  This 
additional  wind  input  contribution  is  parameterized  spectrally  using  the  breaking  probability 
with  wave  scale. 

Once  A(c)  has  been  calculated,  the  peak  breaking  probability  is  obtained  using  equation  (18). 
Fig.  14  above  shows  the  variation  of  spectral  peak  breaking  probability  against  inverse  wave  age 
for  a  range  of  wind  speeds  and  normalised  saturations.  From  this,  it  is  seen  that  the  spectral  peak 
breaking  probability  is  a  strong  function  of  both  wave  age  and  normalized  saturation,  which  is 
expected  based  on  equation  (16)  and  the  relationship  between  normalized  saturation  and  wave 
age.  Further,  there  is  no  indication  of  significant  wind  speed  dependence,  suggesting  that 
breaking  probability  is  much  more  a  function  of  sea  state  than  wind  speed. 

The  self-consistency  of  the  breaking  probability  parameterization  used  as  part  of  the  input  source 
term  is  verified  from  the  computed  breaking  properties  output  from  the  model  for  the  dominant 
(spectral  peak)  waves,  where  such  data  is  available.  Fig.  15  below  shows  close  agreement  for  the 
breaking  probability  against  normalized  spectral  saturation  between  our  parametric  data-fit 
curve,  our  model  forecasts  using  equation  (18),  and  observation  of  breaking  probability. 
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Figure  15.  The  observed  breaking  probabilities  (colored  circles)  are  plotted  against  observed 
normalized  saturation,  as  well  as  our  computed  breaking  probabilities  (open  stars)  from  the  model 
using  equation  (18).  The  model  forecasts  show  a  close  correspondence  with  the  observed  data  and 
the  data-fit  parametric  curve  (solid  line)  also  shown  earlier  in  Fig.  2.  This  validates  the  parametric 
relationship  assumed  in  our  model  to  quantify  the  additional  contribution  of  wave  breaking  to  the 
computed  wind  stress  shown  in  Fig.  10  above. 


Effective  breaking  strength  beff 

As  discussed  earlier  in  section  2.1.4,  we  compiled  a  plot  of  the  variation  of  beff  with  wave  age 
Cp/Uio  for  different  sea  state  conditions.  This  is  shown  in  Fig.  16  below  and  includes  a  best-fit 
curve  for  the  observed  trend.  Also  shown  in  this  figure  are  our  corresponding  model  results 

based  on  our  computed  S^sc(c)  and  A(c)  results,  as  reported  in  Zappa  et  al.  (2016).  The 
corresponding  model  results  of  Sutherland  and  Melville  (2015)  are  also  shown. 


Wave  Age.  Cp  /  u. 


Figure  16.  Variation  of  the  measured  mean  spectral  breaking  strength  b  with  the  wave  age  of  the 

eff 

wind  sea.  The  data  sets  used  are  shown  in  the  right  hand  legend.  The  modeled  results  are  shown  in 
the  left  hand  legend,  including  our  results.  The  solid  black  curve  is  an  optimal  parametric  fit  [from 
Zappa  et  al.,  2016]. 
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4.5  Benchmark  E  -  Overall  consistency  between  wind  input  and  dissipation  rates 


Driving  the  transition  from  growing  to  mature  seas  is  the  source  term  balance,  shown  below  in 
Fig.  17,  for  both  moderate  winds  (6  m/s)  and  strong  (24  m/s).  A  most  interesting  feature  is  that 
even  though  the  wind  input  to  the  dominant  waves  decreases  to  well  below  the  dissipation  rate  as 
the  wave  speed  approaches  the  wind  speed,  the  dominant  wave  saturation  level  (and  steepness) 
remains  sufficiently  large  for  the  breaking  to  occur  through  nonlinear  wave  group  interactions 
(see  Banner  and  Peirson  (2007)). 


_  Source-term  balance  at  U/c  =  1 .5265.  wind  =  6  m/s 
x  10"*  p 


c  Source-term  balance  at  U/c  =  0.71237,  wind  =  6  m/s 
x  10"6  " 


Source-term  balance  at  U/c  =  0.85965.  wind  =  24  m/s 

p 


Figure  17.  Source  term  balance  for  Uio  =  6  m/s  (upper  panels)  and  24  m/s  (lower  panels)  for  young 
(left  panels)  and  maturing  (right  panels)  seas. 


Our  proposed  form  (2)  for  Sds  is  able  to  account  for  breaking  onset  and  loss  rates  during  the 
transition  to  swell.  This  situation  arises  for  maturing  wind  seas,  where  the  wind  input  to  these 
longest  waves  is  greatly  reduced.  At  this  stage  of  evolution  for  the  6  m/s  wind  speed  case,  the 
dominant  balance  is  between  the  dissipation  and  nonlinear  transfer  (see  the  upper  right  panel  in 
Fig.  17  upper  panels),  where  the  dissipation  is  primarily  background  turbulence  as  the  spectral 
peak  saturation  level  falls  below  the  breaking  threshold. 

We  verified  that  during  the  evolution  at  each  wind  speed,  the  integrated  wind  input  and 
dissipation  rates  closely  tracked  each  other,  with  the  dissipation  rate  approaching  the  wind  input 
rate  as  the  seas  age.  The  integrated  nonlinear  transfer  term  was  zero  at  each  timestep.  Given  that 
the  total  wind  input  to  the  waves  is  less  than  the  total  wind  input  to  the  sea  surface,  and  that  the 
drag  coefficients  are  realistic,  the  modeled  dissipation  rates  appear  to  be  very  plausible. 

Our  wind  input  and  dissipation  rate  source  term  behavior  with  wave  age  showed  interesting 
properties  when  reviewed  over  a  wide  range  of  wind  speeds.  These  are  conveniently  plotted  in  a 
normalized  form  in  Fig.  18  using  ul  cp  to  collapse  the  very  large  dynamic  range  with  wind  speed, 
where  cp  is  the  spectral  peak  phase  speed. 
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Figure  18.  Integrated  source  terms  normalized  by  u*2cp  (upper  left)  wind  input  energy  flux 
(upper  right)  dissipation  rate  (lower  left)  local  dissipation  rate  (lower  right)  ratio  of  local  to  total 
dissipation  rate. 


Fig.  18  shows  the  source  terms  and  their  components,  normalised  by  u*2cp.  Even  though  the 
source  terms  vary  over  four  or  more  orders  of  magnitude,  their  normalised  integrals  come  close 
to  falling  onto  single  curves.  The  non-local  components  of  the  dissipation  rate  do  not  normalize 
as  well  as  the  input  or  local  dissipation  rate,  hence  the  total  dissipation  and  the  fraction  of  the 
local-to-total  dissipation  rate  show  a  greater  spread. 

Fig.  19  below  shows  the  behavior  of  the  evolution  of  the  integrated  source  terms  as  the  seas  age 
at  four  different  wind  speeds.  It  is  seen  that  the  integrated  Sin  and  Sds  terms  initially  increase  then 
decrease  towards  large  wave  age.  They  asymptote  towards  each  other  as  the  wind  sea  matures, 
although  the  input  is  always  greater  than  the  dissipation.  This  suggests  that  the  wind  sea  will 
continue  to  grow,  but  at  an  ever-decreasing  rate.  The  integrated  Sni  term  is  zero  for  all  times,  as 
required.  In  addition  to  the  overall  wind  input  momentum  flux,  as  expressed  by  the  drag 
coefficient  (see  Fig.9  in  section  4.3  above),  another  important  validation  check  can  be  made 
using  the  total  wave  energy  dissipation  rate  in  the  water  column,  reflected  in  the  integrated 
dissipation  rate  Sds. 

For  wave  age  around  Cp/Uio  ~  1.2,  the  winds  and  dominant  waves  are  largely  decoupled 
energetically,  and  the  total  dissipation  rate  Sds  in  the  wave  evolution  equation  mirrors  the  total 
dissipation  rate  in  the  wave  boundary  layer,  which  has  been  investigated  observationally  in 
considerable  detail  in  recent  measurements  (Sutherland  and  Melville,  2015,  Fig. 16).  This  has 
allowed  a  valuable  assessment  of  the  validity  of  our  Sds  model  formulation  results,  which 
correspond  closely  to  the  observations  for  the  available  cases. 


29 


Figure  19.  Evolution  of  the  integrated  source  terms  as  the  seas  age,  for  wind  speeds  of  12,  24,  36 
and  48  m/s. 

4.6  Wave  model  performance  -  Case  study  of  very  young  wind  sea  evolution 

We  investigated  the  performance  of  our  refined  source  terms  for  wave  dissipation  and  wind 
input  integrated  over  the  spectrum  against  the  observed  terms  during  the  young  wind  sea  growth 
episode  in  the  Strait  of  Juan  de  Fuca  reported  by  Schwendeman  et  al.  (2014).  The  results  in  Fig. 
20  compare  the  fetch  evolution  of  our  modeled  integrated  Sin  and  Sas  source  terms  with  those 
observed  (left  panel)  and  the  computed  versus  observed  significant  wave  height  Hs  (right  panel). 
The  close  correspondence  of  the  results  is  very  reassuring. 


Significant  Wave  Height  versus  Fetch 


Figure  20.  Comparisons  of  modeled  and  observed  fetch  evolution  of  the  source  terms  for  the  Strait 
of  Juan  de  Fuca  experiment.  Left  panel  shows  the  spectrally-integrated  modeled  wind  input  S[n 
(red),  observed  total  input  (magenta),  computed  total  dissipation  rate  Sds  (blue),  which  is 

dominated  by  the  wave  breaking  component  S^c(cyan),  and  the  observed  dissipation  rate  (black 

dashed  line).  The  right  panel  shows  the  modeled  significant  wave  height  Hs  (blue)  versus  the  two 
available  observations  (cyan,  red). 
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4.7  Summary  of  wave  model  results 

In  this  study,  a  framework  for  incorporating  wave  breaking  predictions  of  breaking  crest  length 
density  and  strength  of  the  dominant  wind  waves  into  sea  state  forecast  models  has  been 
developed.  It  has  been  validated  under  the  assumption  that  it  will  be  implemented  in  conjunction 
with  the  exact  form  of  the  nonlinear  spectral  transfer  term  Sni,  together  with  a  suitable  wind  input 
source  term  formulation  Sin.  Because  of  the  computational  cost,  our  validations  have  been 
limited  to  idealized  cases,  particularly  duration-limited,  unidirectional  (but  possibly  non¬ 
stationary)  winds.  These  cases  provided  a  cost-effective  but  limited  validation  regime  and 
further  validation  is  needed.  This  is  presently  underway,  as  described  in  section  4.10  below. 
Validation  of  breaking  wave  forecasts  at  very  high  wind  speeds  is  clearly  needed,  but  is 
presently  hampered  by  the  lack  of  observational  breaking  data  for  severe  wind  speeds. 

Overall,  the  model  framework  we  have  developed  and  refined  in  this  study  has  performed  very 
well  over  a  particularly  wide  range  of  wind  speeds  according  to  the  five  benchmarks  designed  to 
provide  a  critical  testbed  of  the  model’s  capabilities.  The  outcome  is  that  the  model  is  in 
principle  ready  for  implementation  into  an  operational  prototype,  as  recommended  in  section 
4.10  below. 


4.8  Sea  spray  forecasts 

As  an  additional  output  from  the  wave  model,  we  developed  a  unified  air-sea  flux 
parameterization  that  incorporates  winds,  waves  and  sea  spray  generation. 

We  developed  a  parameterization  structure  for  the  estimation  of  the  fluxes  of  momentum, 
sensible  heat,  and  latent  heat  over  the  ocean.  The  parameterization  is  principally  intended  for  use 
in  coupled  air-ocean-wave  models.  It  is  based  on  a  combination  of  the  NOAA  COARE  bulk 
flux  parameterization  and  the  Fairall-Banner-Morison  (FBM)  sea  spray  parameterization. 
COARE3.5  is  the  latest  version  (Edson  et  al.,  2013).  It  has  been  fitted  to  15,000  hrs  of  direct 
covariance  flux  data  for  wind  speeds  from  0-25  m/s.  The  original  FBM  model  (documented  in 
Fairall  et  al.,  2009)  is  a  scaling  model  where  sea  spray  is  produced  by  breaking  waves.  Droplet 
production  is  driven  by  energy  dissipated  by  wave  breaking.  Ejection  into  the  atmospheric  flow 
is  scaled  by  wind  speed,  slope  of  the  dominant  wave.  Droplet  evaporation  physics  are  used  to 
compute  the  effects  of  the  spray  on  the  total  (direct  and  spray- mediated)  heat  and  moisture 
fluxes.  The  COARE  and  FBM  models  have  been  upgraded  based  on  the  idealized  wave  model 
calculations  using  the  UNSW  wave  model  described  above.  The  wind  stress  in  the  COARE  data 
and  the  wave  model  are  well-described  by  a  Chamock  parameter,  a,  of  the  form  a  =  A(ut/cp)B . 

The  FBM  model  has  been  extended  to  allow  use  of  the  spectral  version  Sds  (k)  of  the  wave 
breaking  dissipation  rate  source  term.  Thus  the  spray  production  by  each  wavenumber  of  the 
breaking  spectrum  is  computed  and  the  total  spray  production  is  the  sum  over  the  spectrum. 
Three  versions  of  the  code  have  been  developed  -  fully  physical,  partially  parameterized  and 
fully  parameterized  -  depending  on  the  data  available  from  the  parent  coupled  model. 

Examples  of  the  model  output  of  heat  fluxes  for  hurricane  surface  layer  conditions  (SST=29°  C, 
air  T=27°  C,  and  RH=90%)  at  24  m/s  and  48  m/s  are  shown  in  Fig.21.  The  use  of  spectrally- 
resolved  droplet  production  has  a  significant  effect  on  droplet  production  in  this  model.  Earlier 
versions  of  the  model  assumed  all  droplets  were  produced  by  the  total  wave  dissipation  from  the 
dominant  wave  breaking.  The  spectral  version  produces  droplets  at  each  wave  phase  speed 
proportional  to  the  dissipation  rate  associated  with  those  waves;  total  production  is  then 
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computed  as  the  integral  over  the  spectrum.  The  result  is  a  relatively  large  increase  in  the 
droplet  mass  flux  from  younger  waves  (see  Fig.21).  We  currently  maintain  codes  of  the  model  at 
three  different  levels  of  computation. 

1.0  Fully  physical:  Here  all  of  the  inputs  to  the  droplet  algorithm  come  completely  from 
elsewhere  (observations,  coupled  atmosphere  wave  model).  The  inputs  are  used  to  compute 
the  drop  spectrum,  the  spectrum  is  used  in  the  thermodynamic  time  scale  code  to  compute 
the  drop  sensible  and  latent  heat  values,  and  the  feedback  code  computes  the  final  values. 

2.0  Partly  parameterized:  Here  fluxes  are  initially  computed  via  a  bulk  flux  algorithm  to  obtain 
u*.  Wave  properties  (Cwave,  hwave,  and  P)  are  computed  using  mathematic  fits  from  a  wave 
model  (e.g.,  the  UMiami,  URI  or  the  UNSW  models).  The  drop  heat  and  moisture  fluxes  are 
explicitly  computed;  feedback  is  computed. 

3.0  Fully  parameterized:  Here  the  basic  bulk  flux  and  wave  property  inputs  are  computed  from 
parameterizations  as  in  case  2.  However,  the  thermodynamics  are  also  calculated  via  a 
parameterization.  Again,  feedback  is  computed.  This  approach  uses  the  least  code  and  is  the 
most  parameterized. 


Figure  21.  Example  model  calculations  of  sensible  and  latent  heat  flux  at  the  surface  using  the 
UNSW  wave  model  output  at  Uio=24  m/s  (left  panel)  and  Uio=48  m/s  (right  panel).  The  red  and 
cyan  lines  are  the  normal  direct  transfers  of  sensible  and  latent  heat.  The  blue  and  green  lines  are 
the  droplet  contributions. 
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Figure  22.  Ratio  of  the  sea  spray  mass  flux  from  our  sea  state  dependent  model  to  that  from  a 
standard  windspeed-only  model  plotted  against  the  inverse  wave  age  (Uio/cp).  Note  how  the  ratio  of 
the  mass  flux  changes  as  the  inverse  wave  age  increases  from  old  seas  (left  side)  where  breaking  at 
the  spectral  peak  is  minimal,  towards  very  young  seas  (right  side)  where  there  is  frequent  breaking 
of  the  spectral  peak  waves.  The  results  shown  are  for  Uio=24  m/s. 


4.9  Full  coupling  of  the  wave  model  to  the  wind  field 

With  the  very  encouraging  validation  to  date  of  our  ‘semi-coupled’  framework,  it  is  feasible  to 
transition  it  to  a  fully  coupled  version.  This  would  require  a  full  atmospheric  model,  using  a 
suitable  wind  input  growth  rate  parameterization  that  ideally  would  depend  on  a  scale-dependent 
wind  velocity,  such  as  U(a/2),  where  A  is  the  water  wavelength.  This  is  presently  ongoing 
through  a  collaborative  effort  with  NCEP  Marine  Modeling  and  Analysis  Branch. 

4.10  Transition  to  operations 

The  breaking  wave  model  framework  developed  in  this  project  has  been  built  around  the  ‘Exact’ 
NL  source  term  for  Sni,  as  we  sought  to  choose  the  most  accurate  source  terms,  albeit  at  the  cost 
of  computational  speed.  However,  the  approximate  form  of  Sni  in  widespread  operational  use, 
known  as  the  Direct  Interaction  Approximation  (DIA),  at  each  timestep  computes  only  one  of 
the  wave-wave  resonance  interaction  terms  of  the  millions  that  are  computed  by  the  ‘Exact’  NL 
source  term.  Needless  to  say,  the  results  can  be  very  different,  and  we  wanted  to  eliminate  this 
from  our  model  development. 

A  much  closer  approximation  known  as  the  Two-Scale  Approximation  (TSA)  has  been 
developed  recently  (Resio  and  Perrie,  2008).  Its  performance  appears  to  match  the  Exact  NL 
very  closely  and  hence  may  be  a  suitable  replacement  for  Exact  NL  in  operations.  This  is  an 
exciting  and  timely  development  that  needs  careful  examination. 

4.11  Highlights  of  the  SPANDEX  study 

SPANDEX,  the  Spray  Production  and  Dynamics  Experiment,  is  a  laboratory  study  conducted  at 
the  UNSW  Water  Research  Laboratory  in  Manly  Vale,  Sydney,  Australia.  The  goals  of 
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SPANDEX  were  to  illuminate  physical  aspects  of  spume  droplet  production  and  dispersion, 
verify  theoretical  simplifications  presently  used  to  estimate  the  source  function  from  ambient 
droplet  concentration  measurements,  and  examine  the  relationship  between  the  implied  source 
strength  and  forcing  parameters  such  as  wind  speed,  surface  turbulent  stress,  and  wave 
properties.  Our  observations  of  droplet  profiles  give  reasonable  confirmation  of  the  basic  power- 
law  profile  relationship  that  is  commonly  used  to  relate  droplet  concentrations  to  the  surface 
source  strength.  This  essentially  confirms  that,  even  in  a  wind  tunnel,  there  is  a  near  balance 
between  droplet  production  and  removal  by  gravitational  settling.  Phase  Doppler  Anemometry 
observations  revealed  significant  mean  horizontal  and  vertical  slip  velocities  that  were  larger 
closer  to  the  surface.  The  magnitude  seems  too  large  to  be  an  acceleration  time  scale  effect,  so 
we  hypothesize  that  the  droplets  tend  to  be  found  in  vertically  moving  air  masses  (e.g.,  updrafts 
departing  from  wave  crests).  Scaling  of  the  droplet  production  surface  source  strength  proved  to 
be  difficult  to  explain  quantitatively.  The  wind  speed  forcing  varied  only  23%  and  the  stress 
increased  by  a  factor  of  2.2.  Yet,  the  source  strength  increased  by  about  a  factor  of  10.  We 
attempted  to  relate  this  to  an  estimate  of  surface  wave  energy  flux  through  calculation  of  the 
standard  deviation  of  small-scale  surface  disturbance  combined  with  wind  forcing.  This  energy 
flux  index  only  increased  52%  with  the  wind  forcing,  so  it  is  not  clear  that  we  have  characterized 
energy  flux  correctly.  Nonetheless,  a  graph  of  spray  mass  surface  flux  versus  surface  disturbance 
energy  is  quasi-linear  with  a  substantial  threshold. 

Laboratory  Study  of  Sea  Spray  from  Breaking  Waves:  Part  I  -  Profiles  of  Droplet  Microphysical 
Properties 

C.W.  Fairall,  C.J.  Zappa,  S.  Brumer,  M.L.  Banner,  R.P.  Morison,  X.  Yan  and  W.L  Peirson. 

This  laboratory  study  was  performed  to  develop  an  accurate  sea  spray  source  function 
parameterization  through  coincident  observations  of  sea  spray  along  with  wave  breaking, 
turbulent  kinetic  energy  dissipation  rate  and  turbulent  fluxes.  It  continued  the  Spray  Production 
and  Dynamics  Experiment  (SPANDEX)  conducted  at  the  UNSW  Water  Research  Laboratory  in 
Manly  Vale  (NSW,  Australia).  SPANDEX  I  is  described  in  Fairall  et  al.  (2009).  SPANDEX  II 
was  carried  out  during  June  2010  and  featured  several  observational  advances  compared  to 
SPANDEX  I:  direct  measurements  of  ocean-side  turbulence  profiles  and  thermal  imaging  of 
individual  spray  drops.  Spray  drops  were  measured  with  an  optical  array  size  spectrometer. 
Properties  of  the  sea  spray  profiles  as  a  function  of  forcing  were  measured  and  compared  to  a 
theoretical  profile  based  on  balancing  turbulent  upward  transport  and  gravitational  settling. 

A  Laboratory  Study  of  Sea  Spray  from  Breaking  Waves.  Part  II  -  Correlations  with  Wind  and 
Wave  Properties.  C.J.  Zappa,  M.L.  Banner,  C.W.  Fairall,  D.  LeBel,  S.  Brumer,  R.P.  Morison, 

X.  Yan  and  W.L.  Peirson. 

Field  measurements  of  spray  flux  show  significant  variability  based  on  wind  speed  or  friction 
velocity.  The  process  of  spray  flux  is  linked  fundamentally  to  wave  breaking,  as  well  as  wind 
speed.  In  SPANDEX  I,  the  laboratory  experiments  used  surrogates  for  dissipation  due  to 
breaking.  In  SPANDEX  II,  we  performed  the  first  study  investigating  spray  flux  dependence  on 
water-side  dissipation  rate,  rather  than  wind  speed.  From  our  laboratory  wind  wave  tank  study 
with  salt  water  (30  psu)  we  found  that  correlating  with  u*3  showed  a  near-linear  dependence. 
Correlating  with  the  integrated  dissipation  rate  also  showed  a  near-linear  dependence.  If  these 
encouraging  initial  results  is  followed  up  in  a  field  study  and  reveal  a  robust  correlation  of  spray 
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flux  with  the  integrated  dissipation  rate,  it  could  be  possible  to  estimate  sea  spray  flux  reliably 
from  the  dissipation  rate  estimated  from  a  wave  model. 

4.12  Benefits  Analysis  Summary 
4.12.1  Technical  Output 

During  the  course  of  this  project,  several  refereed  publications  have  been  published  in  leading 
journals,  with  several  in  preparation,  as  well  as  a  number  of  conference  and  workshop 
contributions. 

Refereed  papers  published 

C.J.  Zappa,  M.L.  Banner,  R.P.  Morison  and  S.E.  Brumer,  2016:  On  the  variation  of  the  effective 
breaking  strength  in  oceanic  sea  states.  J.  Phys.  Oceanogr.  doi:  http://dx.doi.org/10.1175/JPQ- 
D- 15-0227.1  *  "  ~~ 

M.L.  Banner,  X.  Barthelemy,  F.  Fedele,  M.  Allis,  A.  Benetazzo,  F.  Dias  and  W.L.  Peirson,  2014: 
Linking  reduced  breaking  crest  speeds  to  unsteady  nonlinear  water  wave  group  behavior.  Phys. 
Rev.  Lett.  112,  114502. 

M.L.  Banner,  C.  J.  Zappa,  &  J.  Gemmrich,  2014:  A  Note  on  Phillips'  Spectral  Framework  for 
Ocean  Whitecaps,  J.  Phys.  Oceanography,  44,  1727-1734. 

Peirson,  W.L.,  J.W.  Walker  and  M.L.  Banner,  2014:  On  the  microphysical  behaviour  of  wind- 
forced  water  surfaces  and  consequent  re-aeration.  J.  Fluid  Mech.  743,  399-447. 

Tolman,  H.L.,  M.L.  Banner  &  J.M.  Kaihatu,  2013.  The  NOPP  operational  wave  model  improve¬ 
ment  project.  Ocean  Modelling,  70,  2-10.  http://dx.doi.Org/10.1016/i.ocemod.2012.ll.011 

J.R.  Gemmrich,  C.J.  Zappa,  M.L  Banner  and  R.P.  Morison,  2013.  Wave  breaking  in  developing 
and  mature  seas.  J.  Geophys.  Res.  Oceans,  241,  118,  4542-4552,  doi:  10.1002/jgrc.20334. 

Zappa,  C.  J.,  M.  L.  Banner,  H.  Schultz,  J.  R.  Gemmrich,  R.  P.  Morison,  D.  A.  LeBel,  and  T. 
Dickey  (2012):  An  overview  of  sea  state  conditions  and  air-sea  fluxes  during  RaDyO.  J. 
Geophys.  Res.,  117,  C00H19,  doi:  10. 1029/201 1JC007336. 

Banner,  M.L.  and  R.P.  Morison,  2010:  Refined  source  terms  in  wind  wave  models  with  explicit 
wave  breaking  prediction.  Part  I:  Model  framework  and  validation  against  field  data.  Ocean 
Modell.  doi:10.1016/j.ocemod.2010.01.002. 

Refereed  papers  submitted 

X.  Barthelemy,  M.L.  Banner,  W.L.  Peirson,  F.  Dias  and  M.  Allis,  2015:  On  the  local  properties 
of  highly  nonlinear  unsteady  gravity  water  waves.  Part  1.  Slowdown,  kinematics  and  energetics. 
Submitted  to  J.  Fluid  Mech.  (pdf  available  at  arXiv: 1 508.06001 ) 

X.  Barthelemy,  M.L.  Banner,  W.L.  Peirson,  F.  Fedele,  M.  Allis,  F.  Dias,  2015:  On  the  local 
properties  of  highly  nonlinear  unsteady  gravity  water  waves.  Part  2.  Dynamics  and  onset  of 
breaking.  Submitted  to  J.  Fluid  Mech.  (pdf  available  at  arXiv: 1 508.06002  ) 

A.  Saket,  W.L.  Peirson,  M.L.  Banner,  X.  Barthelemy  and  M.  Allis,  2015:  Wave  breaking  onset 
of  two-dimensional  deep-water  wave  groups  in  the  presence  and  absence  of  wind.  Submitted  to 
J.  Fluid  Mech.  (pdf  available  at  arXiv: 1508.07702  ). 

M.L.  Banner  and  R.P.  Morison,  2016:  On  the  upper  ocean  turbulent  dissipation  rate  due  to  very 
short  breaking  wind-waves.  Submitted  to  Ocean  Modelling,  (pdf  available  at  arXiv:  1 602. 0664 9). 
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Refereed  journal  papers  in  preparation,  to  be  submitted 

R.P.  Morison  &  M.L.  Banner,  2016:  Forecasting  wave  breaking  in  ocean  wind  wave  models. 
Part  n.  Model  performance  over  a  wide  wind  speed  range. 

C.W.  Fairall,  M.L.  Banner,  W.L.  Peirson  and  R.P.  Morison,  2016:  Wave  breaking  and 
turbulence  model  of  spume  droplet  production:  a  physically-based  parameterization  of  the  sea- 
spray  large  droplet  source  function. 

C.W.  Fairall,  M.L.  Banner  and  R.P.  Morison,  2016:  A  unified  air-sea  flux  parameterization 
incorporating  winds,  waves  and  sea  spray  generation. 

C.  J.  Zappa,  C.W.  Fairall,  M.L.  Banner,  W.L.  Peirson,  D.  LeBel  and  R.P.  Morison,  2016:  A 
laboratory  study  of  wind-driven  sea  spray  from  breaking  waves. 

Conference  Research  seminars  presented 

R.P.  Morison  and  M.L.  Banner,  Progress  on  Dangerous  Breaking  Wave  Warnings  from  Spectral 
Wave  Forecast  Models.  14th  International  Workshop  on  Wave  Hindcasting  and  Forecasting,  and 
5th  Coastal  Hazards  Symposium,  Key  West,  Florida,  USA,  Nov  8-13,  2015. 

M.L.  Banner,  X.  Barthelemy,  F.  Fedele,  M.  Allis,  A.  Benetazzo,  F.  Dias,  W.L.  Peirson.  Linking 
reduced  breaking  crest  speeds  to  unsteady  nonlinear  water  wave  group  behavior.  Waves  in 
Shallow  Water  Environments  (WISE)  workshop,  8-12  June  2014,  ECMWF,  Reading,  U.K. 

R.P.  Morison  and  M.L.  Banner.  Incorporating  Breaking  Wave  Predictions  in  Spectral  Wave 
Forecast  Models.  13th  International  Workshop  on  Wave  Hindcasting  and  Forecasting  &  4th 
Coastal  Hazards  Symposium,  Oct.  27-  Nov.l  2013,  Banff,  Canada. 

Fairall,  C.W.,  C.J.  Zappa,  S.  Brumer,  M.L.  Banner,  R.P.  Morison,  X.  Yan  and  W.L  Peirson, 

2012:  Laboratory  study  of  sea  spray  from  breaking  waves:  Part  1  -  Profiles  of  droplet 
microphysical  properties.  18th  AMS  Conference  on  Air-Sea  Interaction,  Boston,  8-13  July,  2012 

R.P.  Morison,  M.L.  Banner,  H.  J.  Alves,  P.P.  Sullivan.  Wind  Wave  Model  Performance  in 
Relaxing  Wind  Seas.  12th  International  Workshop  on  Wave  Hindcasting  and  Forecasting  &  3rd 
Coastal  Hazard  Symposium,  Oct.30  -  Nov.  4,  2011,  Kohala  Coast,  Hawaii. 

4.12.2  Strategic  Implications 

Negotiations  are  underway  with  NCEP/Marine  Modeling  and  Analysis  Branch,  NOAA  to 
investigate  transitioning  this  modeling  framework  for  forecasting  wave  breaking  properties  of 
interest  into  WaveWatch  III  for  expanding  the  range  of  forecast  products  to  enhance  public 
warnings  of  dangerous  breaking  wave  conditions. 
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